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The proton momentum distribution, accessible by deep inelastic neutron scattering, is a very sen- 
sitive probe of the potential of mean force experienced by the protons in hydrogen-bonded systems. 
In this work we introduce a novel estimator for the end to end distribution of the Feynman paths, 
i.e. the Fourier transform of the momentum distribution. In this formulation, free particle and envi- 
ronmental contributions factorize. Moreover, the environmental contribution has a natural analogy 
to a free energy surface in statistical mechanics, facilitating the interpretation of experiments. The 
new formulation is not only conceptually but also computationally advantageous. We illustrate the 
method with applications to an empirical water model, ab-initio ice, and one dimensional model 
systems. 

PACS numbers: 05.10.-a, 61.05.F- 



The behavior of protons and more generally of light nu- 
clei in condensed phases is significantly affected by quan- 
tum effects even at ambient temperatures. The isotopic 
effect in water, the ferroelectric behavior of KDP, and 
the formation of high pressure ice phases, are just a few 
of the relevant phenomena where the quantum behav- 
ior of the nuclei plays a role. To address these issues a 
powerful experimental tool, deep inelastic neutron scat- 
tering (DINS) that measures the momentum distribution 
[IH3] has recently been developed. Quantum effects are 
revealed by strong deviations from the classical Maxwell 
distribution. However interpreting DINS experiments is 
difficult and so far has been based on extensive and chal- 
lenging ab initio molecular dynamics simulations [U [5] . 
While these calculations have shown that good agree- 
ment between theory and experiments is possible, a sim- 
pler way of calculating the momentum distribution needs 
to be found and the link between the experimental data 
and the underlying physics made transparent if DINS is 
to become a standard tool. 

In order to understand the source of this computational 
challenge, let us contrast the expression for the momen- 
tum distribution n(p) and that of the partition function 
Z in terms of the density matrix p(r, r') = |e _,3ff | r'). 
The former may be expressed as: 
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where n(x) = ~ J drdr'S (r 
tion function is given by: 
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Z = J dvp{v,v). (2) 

It can be seen that n(p) involves the off-diagonal ma- 
trix elements while Z is determined solely by diagonal 
terms. In a condensed system the potential energy sur- 
face in which the particles move is in a high dimensional 
space and statistical sampling is the only viable computa- 
tional strategy. This is usually done using the Feynman 
path representation. In this representation, n(x) is an 
end to end distribution of a sum over open paths, while 
closed ones determine Z [7] . Sampling is done on the 
closed paths that specify Z and it is challenging from 
these simulations to estimate the open path distribution 
that determines n(p). 

One approach is to artificially open a fraction of the 
paths [8 . In so doing one has to balance two contra- 
dictory requirements. On one hand the number of open 
paths has to be large enough to obtain good statistics for 
n(x), while on the other hand it cannot be too large as 
the sampling will become incorrect. In this work we in- 
troduce a new expression for n(x) which does not require 
opening the paths and compromises neither sampling ac- 
curacy nor statistics. Following a derivation whose detail 
can be found in the supplementary material we find: 
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where no( x ) = e 2l)h ' 2 is the free particle end to end 
distribution. The function y(r) is arbitrary but for the 
boundary condition y(/3h) — y(0) = 1. In practice, the 
optimal choice is to take y = \ — -S^ since it mini- 
mizes the distance between r(r) and the displaced path 
r(r) + y(r)x. Notice that, for simplicity, Eq. ^ refers 
to a single particle subject to the external potential V[r]. 
Generalization to many-body systems is straightforward 
if exchange effects between identical particles can be ne- 
glected. How to include such effects will be discussed in 
a future publication. Eq. (|3| merits further comment. In 
the calculation of the kinetic energy it has been found 
to be extremely useful to use estimators in which the 
free particle contribution has been explicitly accounted 
for [5] . We expect similar computational advantages from 
the explicit separation of no(x). Furthermore it follows 
from Eq. pi) that, having put Z(0) = Z, we can write 
= 5t1y as a ratio between two partition functions. 

n (x) Z(0) r 

To calculate this ratio or its logarithm U(x) = — In 
standard statistical mechanics methods such as free en- 
ergy perturbation |10| or thermodynamic integration [11] 
may be utilized. 

Using free energy perturbation one may compute: 

U( x ) = -ln^e-K/cf^C^W+j/Wxi-yWr)])^ _ ( 4 ) 

where the average is evaluated using the closed path dis- 
tribution Z(0). 

The free energy perturbation method can only be ap- 
plied to systems with weak quantum effects. For sys- 
tems with strong quantum effects the average is diffi- 
cult to converge and instead we use thermodynamic in- 
tegration. In this scheme U (x) is obtained as an integral 
U(x) = L fix' • F(x') over the mean force, 

F(x') = {^ J^dr V r y[r(r) + y(r)x']y(r)^ (5) 

evaluated at the intermediate distributions Z{pd). In 
this case thermodynamic integration requires opening the 
paths, but it does so in a fully controlled way. Besides 
being rigorous our estimator offers several computational 
advantages. In three dimensions the standard approach 
suffers from poor statistics at short distances due to the 
geometrical r 2 factor, and this is not the case here. By 
averaging over all the particles, statistics can be greatly 
improved. The calculation over different particles is in- 
trinsically parallel and the power of modern comput- 
ers optimally harnessed. Furthermore, in crystals where 
anisotropies are relevant, the dependence of n(p) on the 
momentum direction can be easily evaluated. 

We first test our algorithm on a flexible model for wa- 
ter [12]. The simulation box contains 32 water molecules. 
The temperature is set to be 296K. Both protons and 



oxygens are treated by quantum mechanics, and are rep- 
resented by 64 classical beads. The end to end distribu- 
tion is spherically averaged in water. The quantum effect 
for water at room temperature is relatively small [4 J . This 
allows us to use free energy perturbation Q and com- 
pare the results with open path integral simulation [5]. 
In the latter case, in principle one proton path should 
be opened and all other paths should be closed. How- 
ever, the resulting statistics would be poor. In order to 
boost statistics one proton path per water molecule was 
opened, as it was found that this approximation leads 
to a negligible error in the momentum distribution due 
to the relatively weak interaction between protons be- 
longing to different water molecules [5]. The closed path 
formulation allows one to compute the end to end distri- 
bution without opening any proton path, and therefore 
all the protons can be included in the calculation of the 
end to end distribution without any approximation. We 
show the end to end distribution calculated both from a 
268 ps open path simulation and from a 12 ps closed path 
simulation that utilizes the estimator given by Eq. Q in 
Fig. [T] (a), and the comparison of the potential of mean 
force in Fig. [I] (b). In both simulations, the time step is 
0.24 fs. Two consecutive steps contain highly correlated 
information, and the free energy perturbation estimator 
may be computed every 20 steps. Thus with only a small 
increase in computational overhead in comparison to an 
open path simulation of the same length, the displaced 
path formulation has a large gain in terms of sampling 
this property efficiently. 
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FIG. 1: (color online) Comparison of (a) the end to end distri- 
bution and (b) the potential of mean force in SPC/F2 water. 
In both figures, the red line is computed by a 268ps open path 
integral simulation. The thick blue line is calculated using the 
displaced path estimator Q, with the thickness indicating the 
95% confidence interval. The noise near r — in both insets 
for open path simulation is due to the r 2 weight in the spher- 
ical integration, while the displaced path gives correct small 
r behavior by definition. 

The thermodynamic integration approach given in 
Eq. ([5]) is not only computationally advantageous but 
also provides one with the potential of mean force U(x), 
and its gradient F(x) which are key quantities for inter- 
preting the physics underlying n(p). We first note that 



3 



the kinetic energy K is given by K 



x=0 



• F(x) 

= Ky + Since 3/2/3 is the free particle contribu- 
tion, the non-classical contribution is completely included 
in the excess kinetic energy term Ky, and is determined 
by the zero point curvature of U (x) . Secondly, if the mo- 
mentum distribution of an individual particle is accessi- 
ble (as is possible e.g. in simulations) and the underlying 
potential energy surface is harmonic, the end to end dis- 
tribution follows a Gaussian distribution and the mean 
force is given by a straight line. Any deviation of q • F(x) 
from linearity signals anharmonic behavior along the q 
direction. 

In experiments, the spherically averaged momentum 
distribution is accessible in liquids, and amorphous and 
polycrystalline solids, while the directional distribution is 
also accessible in mono crystalline materials. The latter 
distribution provides more information about the under- 
lying potential energy surface. However, in single crystals 
the total momentum distribution is the sum of the con- 
tributions of individual particles participating in bonds 
with different orientations. As a consequence the differ- 
ence between directional and spherical momentum distri- 
bution is usually very small as shown in the top panel of 
Fig. [2] This figure is based on an anisotropic harmonic 
model [13] with three distinct principal frequencies that 
is fit to the ab initio path integral data for ice Ih [?]. 
The bottom panel of the same figure clearly shows that 
the distinction between the spherical and directional dis- 
tributions is enhanced when comparing the mean forces. 
It is therefore of great interest to link directly the mean 
force to the experimental data, i.e. to the Compton pro- 
file J(q, y) = J n(p)5(y — p • q)rfp where q indicates the 
direction of the neutron detector [2|. One finds with a 
derivation provided in the supplemental material that the 
mean force is related to the Compton profile by: 
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In the bottom panel of Fig. [2]the slope of the mean force, 
either spherical or directional, at r = is equal to the 
excess kinetic energy Ky divided by the constant 
This is an exact result that originates from the symme- 
try property of ice Ih. In general the spherical and direc- 
tional mean force can have different slopes at r = 0. The 
deviation of the spherical and directional forces from lin- 
earity at finite r results from the averaging process and is 
not a sign of anharmonicity. Thus in the interpretation 
of the experimental Compton profile, which results from 
the contribution of many particles, one must distinguish 
the case of an anisotropic harmonic potential energy sur- 
face from that of an anharmonic potential energy surface. 
To the best of our knowledge the procedure that is cur- 
rently adopted to fit the experimental data |31 [H] does 
not separate well anisotropic and anharmonic effects. We 
propose here an alternative approach in which the mean 



force is associated to the experimental Compton profile 
according to Eq. The projections of the mean force 
along different directions are then fitted to an anisotropic 
harmonic model averaged as required by the crystal sym- 
metry. Any systematic deviation from experiment of the 
mean force originating from the harmonic contribution, 
can then be associated to anharmonicity and used to fur- 
ther refine the underlying model potential energy surface. 





FIG. 2: (color online) Top panel: the momentum distribu- 
tion of the protons in ice Ih resulting from an anisotropic 
harmonic model (see text). Both the spherical and the direc- 
tional distribution along the c-axis are shown. Bottom panel: 
the corresponding spherical and directional mean force pro- 
jected along the c-axis. The curves are plotted as a function 
of the end to end distance. The mean force enhances the 
differences between spherical and directional distributions. 

The framework introduced here may be also utilized to 
provide insight to the investigation of anharmonic sys- 
tems. Consider for example a particle with the pro- 
ton mass subject to a model double well ID-potential. 
V = =f^z 2 + Aexp(-^) with w = 1578K, and £ = 

0.094A. A characterizes the barrier height and is set to 
be 1263K, 3789K, and 6315K, respectively. These param- 
eters mimic different tunneling regimes for protons along 
a hydrogen bond [5J [TS] . The temperature is set to be 
30K. At this temperature the behavior of the systems 
is dominated by the ground-state, and the end to end 
distribution can be approximated by the overlap integral 
n(x) = J dzip(z)ip(z + x) where ip(z) is the ground-state 
wavefunction and F(x) = — ^ lnn(x). In Fig. [3] we can 
see how qualitatively different the mean force can be in 
the three cases. One goes from a fully monotonic behav- 
ior for A — 1263K which is a model for a low energy 
barrier hydrogen bond [16] , to the strongly non mono- 
tonic mean forces for A — 3789K, A = 6315K where the 
tunneling states lie below the barrier height. Addition- 
ally, it is not very difficult to relate features of the mean 
force to the underlying effective potential. 

It is also instructive to study F(x) as a function of 
temperature when the higher states are mixed in the 
density matrix. This is done in Fig. [4] for the double 
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FIG. 3: (color online) (a) The mean force corresponding to 
a double well model at T = 30K, for different barrier heights 
A = 1263K (black solid line), A = 3789K (red dashed line), 
and A — 6315K (blue dot-dashed line), (b) Potential energy 
surface for A — 1263K (blue solid line), and the first five 
energy levels (red dashed line), (c) (d) the same as (b), but 
with A = 3789K and A = 6315K respectively. 
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FIG. 4: (color online) The mean force corresponding to a 
double well model at A = 3789K for different temperatures 
100K (red solid line), 300K (blue triangle), 500K (black dot- 
dashed line), 1000K (magenta dashed line), and 2000K (blue 
cross). 



well potential with A = 3789K. For temperatures in the 
100 — 500K range, the behavior is dominated by the two 
lowest eigenstates. The slope of F(x) at small x, which is 
proportional to the excess kinetic energy Ky , shows little 



dependence on T. It can be shown with detailed analysis 
that this is a generic feature of two level tunneling sys- 
tems. Other characters seen in Fig. [3] in the same range 
of temperatures, such as the more pronounced kink at 
intermediate x and the enhanced softening of the mean 
force at large x, derive from the odd symmetry of the 
first excited state contribution. Eventually at higher T 
the kink in F(x) disappears as the mean force progres- 
sively resumes linear behavior with a slope that tends to 
zero as high temperature classical limit is reached. 

In this work, we develop a novel displaced path for- 
malism for the calculation of momentum distribution of 
quantum particles. The algorithm is rigorous and com- 
putationally advantageous. The new formulation intro- 
duces in a natural way a potential of mean force which 
is a quantity that very clearly illuminates the physics be- 
hind n(p) and can be used to further understand and 
analyze experimental and theoretical results. 

This work is partially supported by NSF under grant 
CHE-0956500 and by DOE under grant DE-FG02- 
05ER46201 (LL and RC). 



* Present address: Department of Chemistry, Columbia 
University, New Yo rk NY 10027 



[9 

[10 

[H 
[12 

[13 

[14 

[15 
[16 



' Electronic address: rcar@princeton.edu 



G.F. Reiter, J. Mayers, and P. Platzman, Phys. Rev. 
Lett. 89, 135505 (2002). 

G.F. Reiter, J. Li, J. Mayers, T. Abdul-Redah, and 
P. Platzman, Braz. J. Phys. 34, 142 (2004). 

C. Andreani, D. Colognesi, J. Mayers, G.F. Reiter, and 
R. Senesi, Advances in Physics 54, 377 (2005). 

J. A. Morrone and R. Car, Phys. Rev. Lett. 101, 017801 
(2008). 

J. A. Morrone, L. Lin, and R. Car, J. Chem. Phys. 130, 
204511 (2009). 

D. Ceperley, Reviews of Modern Physics 67, 279 (1995). 
D. Ceperley and E. Pollock, Can. J. Phys. 65, 1416 
(1987). 

J. A. Morrone, V. Srinivasan, D. Sebastiani, and R. Car, 
J. Chem. Phys. 126, 234504 (2007). 

M. Herman, E. Bruskin, and B. Berne, J. Chem. Phys. 
76, 5150 (1982). 

R. Zwanzig, J. Chem. Phys. 22 (1954). 

J. Kirkwood, J. Chem. Phys. 3, 300 (1935). 

J. Lobaugh and G. A. Voth, J. Chem. Phys. 106, 2400 

(1997) . 

L. Lin, J. A. Morrone, R. Car, and M. Parrinello (2010), 
in preparation. 

A. Pietropaolo, R. Senesi, C. Andreani, A. Botti, 
M.A. Ricci, and F. Bruni, Phys. Rev. Lett. 100, 127802 
(2008). 

M. Benoit, D. Marx, and M. Parrinello, Nature 392, 258 

(1998) . 

M. Benoit and D. Marx, ChemPhysChem 6, 1738 (2005). 



